AD-A150  792  LOCAL  COVARIANCE  FUNCTIONS  AND  DENSITY  DISTRIBUTIONS 
(U>  OHIO  STATE  UNIV  RESEARCH  FOUNDATION  COLUMBUS 
R  FORSBERG  JUN  84  OSU/DGSS-356  AFGL-TR-84-0214 
UNCLASSIFIED  F19628-82-K-0022  F/G  8/5 


AD-A150  792 


• 


AFGL-TR- 84-0214 


> 


LOCAL  COVARIANCE  FUNCTIONS  AND  DENSITY  DISTRIBUTIONS 


* 


Rene  Forsberg 


The  Ohio  State  University  __ 

Research  Foundation  [J 

Columbus,  Ohio  43212 

r. 

£ 


dune,  1984 


Scientific  Report  No.  6 

Approved  for  public  release;  distribution  unlimited 


AIR  FORCE  GEOPHYSICS  LABORATORY 
AIR  FORCE  SYSTEMS  COMMAND 
UNITED  STATES  AIR  FORCE 
HANSCOM  AFB,  MASSACHUSETTS  01731 


CONTRACTOR  REPORTS 


This  technical  report  has  been  reviewed  and  Is  approved  for  publication. 


CHRISTOPHER  JEKw 
Contract  Manager 


THOMAS  P.  ROONEY  7 

Chief,  Geodesy  &  Gravity  Branch 


FOR  THE  COMMANDER 


DONALD  H.  ECXHARDT 
Director 

Earth  Sciences  Division 


This  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  Is 
releasable  to  the  National  Technical  Information  Service  (NTIS) . 


Qualified  requesters  may  obtain  additional  copies  from  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical 
Information  Service. 


If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing 
list,  or  If  the  addressee  Is  no  longer  employed  by  your  organization,  please 
notify  AFGL/DAA,  Hanscora  AFB,  MA  01731.  This  will  assist  us  In  maintaining 
a  current  mailing  list. 


Unclassified _ 

SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Wlitn  Data  Entarad) 


REPORT  DOCUMENTATION  PAGE 


It.  REPORT  NUMBER 


AFGL-TR-84-0214 


1 4.  TITLE  (and  Sublllla) 


3  A/'P  READ  INSTRUCTIONS 

_ BEFORE  COMPLETING  FORM 

2.  GOVT  ACCESSION  NO.  3  RECIPIENT'S  CATALOG  NUMBER 


5.  TYPE  OF  REPORT  A  PERIOO  COVERED 


Local  Covariance  Functions  and 
Density  Distributions 


17.  AUTHORC*; 


Rene  Forsberg 


Scientific  Report  No.  6 

S.  PERFORMING  ORG.  REPORT  NUMBER 

OSU/DGSS-356 _ 

8.  CONTRACT  OR  GRANT  NUMBERS) 

F19628-82-K-0022 


"•  • 


•  • 


Ji  -  -  — 1  s  —a— — V  i -i  '4 

•  • 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS  10.  PROGRAM  ELEMENT.  PROJECT.  T ASK 

T,  ~  .  .  ...  .  _  AREA  &  WORK  UNIT  NUMBERS 

The  Ohio  State  University  _ 

Research  Foundation  oUorinr 

1958  Neil  Avenue,  Columbus,  Ohio  43210 

II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  '2-  REPORT  DATE 

Air  Force  Geophysics  Laboratory  June  1984 

Hanscom  AFB,  Massachusetts  01731  o.  number  of  pages 

Monitor/Christopher  Jekeli/LWG  56 

14.  MONITORING  AGENCY  NAME  A  AOORESSflf  diltarant  from  Controlling  Oltica)  IS.  SECURITY  CLASS,  (ol  thia  raport) 

Unclassified 


<5*.  DECLASSIFICATION/  DOWNGRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (ol  thia  Raport) 


Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (of  thm  abatract  antarmd  in  Block  20,  it  different  trom  Report) 


ti.  SUPPLEMENTARY  notes 


fir  K.EY  WOROS  (Contlnum  on  reveree  aida  ii  neceaaary  and  identify  by  bloc*  number) 

*t>su  t**x<*/* -t~  .■ 

Gravity,  Covariance  Functions,  Density  Distributions,  Collocation, 


20.  ABSTRACT  ( Continue  on  raver  ae  aida  It  nacaaaary  and  idmntlfY  by  block  number) 

>  The  relationship  between  covariance  functions  and  the  density  anomaly  dis¬ 
tributions  generating  the  gravity  field  is  studied  using  the  ensemble  averaging 
theorem,  yielding  interpretations  of  common  well-known  covariance  functions  in 
terms  of  simplified  statistical  mass  models. 

Emphasis  is  on  application  for  local  gravity  field  modelling  (e.g. 
assuming  a  high  degree-and-order  spherical  harmonic  reference  field  to  be 
used),  within  trie  framework  of  the  planar  approximation.  The  very  simple  rela- 


DD  /A 


1473  EDITION  OF  I  NOV  AS  IS  OBSOLETE 


■'«  * 

*«*  * 
.*  V  V  \ 


. 

‘VVVr'A*  V*' 


VA--V-V-V.  > 
'«*-•••  *.  ••  * 


Unclassified _ 

S E CuRtTY  CLASSIFICATION  OF  THIS  PAGE  '(When  Data  Entered) 


•  '  S*  *.W  Va' 


••-V.  *\  ■* 


.>  .■  .■ 


SECURITY  CLASSIFICAIIUN  OF  THIS  PAOEflWnn  D*t*  Entmd)  . 


^tionship  existing  between  the  planar  power  spectrum  and  the  degree  variances 
are  treated  in  detail,  and  it  is  outlined  how  the  shape  of  the  power  spectrum 
may  be  used  as  a  geophysical  inversion  tool,  to  yield  depths  to  density  con¬ 
trast  interfaces  within  the  earth. 

As  a  special  application  of  the  simple  covariance  functions  associated 
with  statistical  mass  distributions,  it  is  shown  how  least-squares  collocation 
may  be  interpreted  as  generalized  point  mass  modelling. 

Finally,  formulas  are  given  for  practically  applicable  local  multi-layer 
covariance  models  ('Compensated  Poisson  model *77  and  gravity  anomalies  in  a 
number  of  sample  areas  in  the  United  States  are  analyzed  to  yield  empirical  co- 
variance  functions,  power  spectra  and  degree-variances. \  Multi-layer  covariance 
models  are  given  for  each  area,  yielding  quite J  similar  results,  and 
demonstrating  the  generally  logarithmic  decay  (-f3)  of  the  degree  variances, 
even  for  degrees  up  to  *=5400.  \ls 


UNCLASSIFIED  _ 

SECURITY  CLASSIFICATIOM  TkJ"  P*GE(Whm  DU*  Enl*r*d) 


This  report  was  prepared  by  Rene  Forsberg,  Geodetic  Institute,  Denmark, 
and  Research  Associate,  Department  of  Geodetic  Science  and  Surveying,  The  Ohio 
State  University,  under  Air  Force  Contract  No.  F19628-82-K-0022,  The  Ohio  State 
University  Project  NO.  714274.  The  contract  covering  this  research  Is  administered 
by  the  Air  Force  Geophysics  Laboratory,  Hanscom  Air  Force  Base,  Massachusetts, 
with  Dr.  Christopher  Jekeli,  Scientific  Program  Officer. 

Certain  computer  funds  used  in  this  study  were  supplied  by  the  Instruction 
and  Research  Computer  Center  through  the  Department  of  Geodetic  Science  and  Surveying 

This  report  is  a  contribution  to  a  gravity  field  modelling  project  supported 
by  NATO  grant  320/82.  Travel  support  for  the  author's  stay  at  The  Ohio  State 
University  have  been  provided  by  a  NATO  Science  Fellowship  Grant. 

The  author  wishes  to  thank  Dr.  Richard  H.  Rapp  for  hosting  me  at  OSU  and  for 
many  fruitful  discussions,  Laura  Brumfield  for  typing  the  report,  and  finally  a 
thanks  to  the  many  graduate  students  of  the  Department,  who  have  assisted  me  through 
valuable  discussions  and  practical  advice. 


Accession  For 

NTIS  GRAJfcl 

US 

DTIC  TAB 

fS 

Unannounced 

□ 

Justification _ 

By_ 


I  Distribution/ 


Dist 


Availability  Codes 
Avail  and/or 
Special 


I 


-iv- 

Con tents 


1.  Introduction  1 

2.  Spectral  Analysis  of  Covariance  Functions  3 

2.1  The  Planar  Approximation  4 

2.2  Hankel  Transforms  6 

2.3  Covariance  Function  and  Power  Spectrum  7 

2.4  Relationship  Between  the  Degree- Variances  and  the  Power  8 

Spectrum 

3.  The  Relationship  Between  Covariance  Functions  and  Mass  10 

Distributions 

3.1  Deterministic  Power  Spectra  of  Elementary  Masses  10 

3.2  Statistical  Mass  Distributions  13 

3.3  Relationship  to  Traditional  Covariance  Functions  17 

4.  Collocation  as  Generalized  Point  Mass  Modelling  22 

4.1  Covariances  as  Gravimetric  Effects  22 

4.2  Application  to  Least  Squares  Collocation  24 

5.  Application  of  Simple  Covariance  Models  for  Local  Covariance  27 

Functions 

5.1  Use  of  Multi-Layer  Models  27 

5.2  Modifications  to  Eliminate  the  Influence  of  Low  Frequencies  29 

6.  Examples  of  Local  Empirical  Covariance  Functions  and  Power  Spectra  35 

7.  Summary  and  Conclusions  48 

References  50 


v; 


1.  Introduction 

Covariance  functions  for  the  gravity  field  form  the  basis  for  optimal 
gravity  field  approximation  techniques  such  as  collocation  and  Wiener  filtering, 
andarealso  indispensable  for  error  studies  and  similar  topics  relating  to  inte¬ 
gral  formulas  etc.  Considering  collocation,  it  is  generally  stated  that  the 
result  of  application  of  collocation  techniques  is  insensitive  to  the  actual 
choice  of  covariance  function  parameters,  whereas  the  error  estimates  are  critically 
dependent  on  these  parameters.  This  is  a  somewhat  dangerous  statement:  experience 
has  shown  it  indeed  to  be  the  case  for  e.g.  gravity  interpolation,  but  when 
predicting  a  different  "order"  of  gravity  field  quantities,  e.g.  geoid  undulations 
from  gravity,  it  is  definitely  not  the  case:  the  predicted  geoid  may  be  too 
rough  or  too  smooth  depending  on  whether  the  implied  covariance  model  has  too 
much  or  too  little  power  in  the  longer  wavelengths. 

There  is  therefore  a  definite  need  for  "good"  covariance  models,  i.e.  parametrical 
models  giving  a  good  fit  to  empirical  data  from  a  given  area.  Many  "simple" 
models  for  parametric  covariance  functions  have  been  suggested  and  evaluated 
over  the  years;  for  a  review  see  e.g.  Moritz  (1980).  Possible  misfits  of  the 
"simple"  models  may  be  reduced  by  "multiple"  models,  as  e.g.  suggested  by  Jekeli 
(1978),  Jordan  (1978)  and  many  others.  These  simple  and  improved  covariance 
models  may  be  interpreted  in  terms  of  statistical  properties  of  the  density 
anomal ies  generating  the  anomalous  gravity  field.  In  the  author's  opinion,  such 
interpretation  models,  outlined  in  the  sequel,  will  principally  be  of  value 
in  clarifying  the  relationship  between  the  various  types  of  simple  covariance 
models  in  use,  and  for  providing  geophysically  reasonable  improved  covariance 
function  models. 

Another  -  and  very  important  application  -  is  the  use  of  "covariance  inter¬ 
pretation"  as  a  geophysical  exploration  method,  to  obtain  depths  to  density  anom¬ 
alies.  Such  statistical  methods  have  been  dominating  in  potential  field  geophysics 
for  the  last  decade  or  so,  and  have  been  applied  successfully  and  routinely 
for  aeromagnetic  data.  For  gravity,  however,  the  successes  have  been  more  limited, 
reflecting  partly  the  usually  much  less  dense  gravity  coverage  available  and 
-  especially  -  that  the  sources  of  the  gravity  field  variations  are  more  "spread 
out"  and  less  applicable  for  the  "white  noise  layer"  descriptions  which  magnetic 
sources  often  seem  to  follow,  e.g.  when  unmagnetic  sediments  are  overlying  a 
crystalline  basement. 
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In  the  sequel,  the  relationship  between  “statistical"  density  anomaly  dis¬ 
tributions  and  associated  covariance  functions  will  be  outlined.  As  a  supplemental 
remark,  the  formulation  expressing  collocation  as  simple  generalized  point  mass  model¬ 
ing  will  be  mentioned,  pointing  out  the  simple  relationship  existing  between  collo¬ 
cation  and  geophysical  inversion  techniques.  Through  the  use  of  such  geophysical- 
oriented  approaches  I  think  that  some  of  the  basics  of  collocation  are  "demystified", 
giving  e.g.  physical  significance  to  such  things  as  the  depth  to  the  Bjerhammar 
sphere. 

In  this  presentation,  I  will  focus  on  local  covariance  functions,  applicable 
only  for  a  given  area  or  geologic  province.  Certainly  the  earth's  gravity  field 
is  in  no  way  a  stationary  process,  and  local  empirical  covariance  must  be  used 
as  a  guideline  for  choosing  optimum  parametric  covariance  models  (or,  rather, 
the  “best”  Hilbert  space  kernels)  for  such  an  area.  For  local  applications,  global 
information  available  through  spherical  harmonic  expansions  should  always  be  utilized 
by  subtracting  such  a  reference  field.  The  remaining  "residual"  field  will  have 
less  variance  and  shorter  correlation  length,  and  especially  when  using  the  available 
high  degree-and-order  expansions  complete  to  degree  180 (e.g.  "Rapp-180",  Rapp 
1981)  the  "variations"  will  be  so  local,  that  the  flat-earth  approximation  becomes 
completely  justifiable. 

This  planar  approximation  will  be  the  object  for  the  present  work.  To  exag¬ 
gerate:  by  the  excellent  spherical  harmonic  models  now  available,  only  local 
problems  remain  in  physical  geodesy!  And  locally  everybody  knows  the  earth  is  flat. 

In  the  planar  approximation,  the  discrete  spherical  harmonic  spectrum  will  be 
replaced  by  the  continuous  Fourier  spectrum.  Although  no  simple  transition  exists 
between  these  spectra,  a  remarkably  simple  and  accurate  transition  exists  between 
the  power  spectrum  and  the  degree-variances.  This  simple  relationship  will  be 
treated  in  some  detail  in  the  next  section.  Due  to  this  simple  relationship,  re¬ 
sults  for  the  planar  approximation  given  in  the  sequel  are  more  or  less  directly 
transferrable  to  a  spherical  earth. 

At  the  end  of  the  report,  some  results  for  actual  data  in  the  U.S.  will  be 
used  as  illustrations,  to  show  some  of  the  applications  and  limitations  of  the 
simple  covariance  models  and  implied  density  anomaly  distributions, and  to  compare 
the  "local"  results  to  conventional  global  covariance  functions. 


2.  Spectral  Analysis  of  Covariance  Functions 

The  anomalous  potential  of  the  earth  defined  as  the  difference  between  the 
actual  geopotential  Wand  a  normal  “ellipsoidal"  potential  U,  may  be  expanded 
in  fully  normalized  spherical  harmonics  as 

"P  om(  si  n  $)  cos  mx 

“ini  tm  ' 

T<r’ *' ‘i  ‘  (m-0)  <2 

.  P__( si n  <(>)  sin  mx 
W  (m  <  0) 

The  function  T  will  be  harmonic  (v2T=0)  outside  a  sphere  of  radius  R,  usually 
taken  as  a  mean  earth  radius  in  available  solutions.  The  spherical  harmonic 
"power  spectrum"  is  the  well-known  degree-variances 
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If  the  covariance  function  K  of  T  is  assumed  to  be  isotropic ,K(P,Q)=K(<|i ,rp,  r^) 
we  have  the  following  well-known  expression  for  the  spatial  covariances 


K(P,  Q)  -  l  o 
1=2  4 


’jlL|*+ 1 

rPrQ 


Pfc(cos  ’ll) 


see  e.g.  Moritz  (1980).  On  the  reference  sphere,  K  and  are  thus  related  by 
a  Legendre  transform,  and 

„  =  2l+l  rn  ut,\  n  ,i.  \  10 


J  K(o0  P4(cos  ’ll)  sin«pdtp 


Figure  1 


The  terms  K  and  o  relate  to  the  anomalous  potential  T>  and  thus  by  division 

Xf 

with  normal  gravity  y  to  the  height  anomaly  (spatial  geold  undulation)  c 


Z  (r.  *  ,  X  )  *  - 


(2.5) 


For  gravity  anomalies 


(2.6) 


similarly  C  and  c^  are  traditionally  used  for  covariance  function  and  degree- 
variances  respectively. 

C  may  be  expressed  by  K  using  covariance  propagation: 


c<p.  «>  *  CCk(’-  •> 


(2.7) 


where  the  dots  indicate  that  the  gravity  anomaly  operators  LAg  should  be  applied 

each  to  one  of  the  variables  in  K(P,  Q).  In  the  spectral  domain,  the  eigenvalue 

£-1 

of  operator  (2.6)  is  well  known  to  be  (Heiskanen  and  Moritz,  1967,  p.  97), 
i.e.  for  the  spherical  harmonic  coefficients 


(2.8) 


and  thus 
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(2.9) 


2.1  The  Planar  Approximation 

The  previous  topics  should  be  familiar  to  most  geodesists.  Less  familiar 
is  often  the  use  of  Fourier  transformations.  The  Fourier  techniques  have  as 
advantages  their  mathematical  simplicity,  and  for  practical  applications,  the 
efficient  transformations  algorithms  available  thorugh  FFT  -  Fast  Fourier  Transform. 


The  drawback  is  that  they  can  only  be  used  local 1y,  since  the  flat-earth  approxi¬ 
mation  forms  the  basis.  However,  with  the  available  high-degree  and  order  spherical 
harmonic  expansions  of  the  geopotential  complete  to  degree  180,  (e.g.  Rapp,  1981), 
the  residual  field 


T'  =  T  -  T 


180X180 


(2.10) 


may  be  very  suitable  for  flat-earth  methods:  T*  will  have  a  covariance  function 
of  short  correlation  length  (typically  30-50  km),  much  smaller  than  distances 
in  which  earth  curvature  effects  need  be  taken  into  account. 

In  the  flat-earth  planar  approximation  the  reference  sphere  will  be  replaced 
by  a  tangential  plane  n.  An  XYZ-system  will  be  used  in  n,  with  X  positive  east, 

Y  positive  north  and  Z  positive  upwards  (Figure  1).  The  Fourier  transformation 
pair  is  defined  through 


1(U,  v,  z)  •  J  Tlx,  y,  z)e"*(ux+vy^dxdy 

n 

T(x,  y,  z)  =  /  T(u,  v,  zje1 ^UX+vy^dudv 


(2.11) 


where  the  spatial  frequencies  or  wave  numbers  will  be  termed  (u,  v)  in  this  report, 
in  accordance  with  the  notation  of  Papoulis  (1968).  The  spectra  at  various  alti¬ 
tudes  z  are  related  through 


T(u,  v,  z)  =  T(u,  v)e 


■cuZ 


*  /u2+v2 


(2.12) 


where  T(u,  v)  is  the  spectrum  at  ihe  reference  plane  n  and  u,  the  so-called  radial 
wave  number. 

Spectra  of  other  gravity  field  quantities  are  easily  found  from  their  "defining" 
operators  e.g. 


°Gravity  anomalies  Ag  = 


.  a  i 


(2.13) 


°N-S  deflections 


°E-W  deflections 


_ iT' 

3r  r 


Ag  =  («w  -  -^)T  =  -wT 
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n  = 


yR  3<P  ‘ 

l  -  T 

Y 

(2.14) 

.  1 

11. 

;  =  -iur 

(2.15) 

Y  COS(j)  R 

3  A  * 

n  Y  T 

The  second  term  in  (2.13)  -  the  "indirect"  effect  -  may  often  be  neglected  when 
using  a  180x180  reference  field:  typical  errors  for  the  reference  field  are  at 
most  a  few  meters  in  e»  corresponding  to  a  residual  indirect  effect  of  a  fraction 
of  a  mgal .  In  this  case,  an  important  relationship  between  the  variances  of  de¬ 
flections  and  gravity  may  be  derived.  For  an  arbitrary  function  f  Parsevals  equation 
says 


j  f 2 ( x ,  y)  dxdy  =  f  f2(u,  v)  dudv 


(2.16) 


In  terms  of  variance,  we  have  for  centered  quantities 


2  =  E{ Ag2 }  = -^EUg2}  =  ^  u,2  E{T2} 


Ag 


4  tt  ' 


(2.17) 


1  ,u2  v2 


1  u,2 


2  +  a2  =  -7T-1  (rr  +  -V)  E{T2}  =  ST7  E{T2}  =  -4  a 

£  n  4tt  2  V2  Y2'  4irz  y2  J  x2  Ag 


(2.18) 


For  an  isotropic  field  a2=a2,  and  thus  o2=4- t  a2  ,  in  terms  of  r.m.s.  variation 

s  n  ?  2r  ig 

corresponding  to  c.  6.7  mgal/arcsec. 


2.2  Hankel  Transforms 

When  Fourier  analyzing  radially  symmetric  functions,  so-called  Hankel  transforms 
are  obtained. 

Let  f(s)  =  f(x,  y),  s  *  /x^+y7  be  a  radially  symmetric  function  in  n.  Then 
the  Hankel  transform  pair  is  (Papoulis,  1968): 


(2.19) 


fU)  =  /*  sf (s)J0(ujs)ds 
0 

f(s)  =  u>  ?■  (a> )  J  (uis)dw  (2.20) 

0  u 

where  J0(-)  is  the  Bessel  function  of  order  zero,  the  Hankel  transform  and  the 
Fourier  transform  of  f  being  essentially  equivalent  since 


f(u,  v)  =  27rf(ii)),  in  =  /u^+v2 


(2.21) 


To  give  two  important  examples  of  Hankel  transforms,  consider  the  reciprocal  dis¬ 
tance  function  ~  =  (x2+ y2+z2)-'5.  The  transform  pairs 


^  1_  e-u)Z 

r  oi 


(2.22) 


1  -a)Z 

—  e 


(2.23) 


may  e.g.  be  found  in  (Papoulis,  1968,  p.  145). 

For  practical  computations,  efficient  numerical  algorithms  ("Fast  Hankel 
Transform" ) exist  ,  building  on  completely  different  principles  than  the  FFT- 
algorithm  (Johansen  and  Sdrensen,  1979). 


2.3  Covariance  Function  and  Power  Spectrum 

In  the  planar  approximation  the  spatial  covariance  function  K  becomes  a  function 
defined  in  the  space  above  the  reference  plane  n: 


rp,  rq)  =  K(s,  zx,  z2)  (2.24) 

where  ss0R  is  the  distance  in  the  plane.  Consider  for  a  moment  the  covariance 
function  at  the  reference  level  rp=rq  =  R,  Zp=ZQs0.  The  discrete  degree-variance 
spectrum  (2.2)  will  in  the  planar  case  be  replaced  by  a  continuous  function  - 
the  power  spectrum  (or  rather  power  spectral  density  -  psd) 
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^_(u,  v)  ■  E(T(u,  v)2}  3  K(u,  v) 


(2.25) 


The  power  spectrum  and  covariance  function  are  thus  related  through  a  Hankel  trans¬ 
form  when  the  covariance  function  Is  assumed  to  be  Isotropic 


4>tt(w)  s  2itK(s)  3  2trK(i^R) 


(2.26) 


Spatial  covariance  expressions,  corresponding  to  (2.3),  may  be  obtained  using 


the  upward  continuation  operator  e”“z  twice: 


K(s,  zp,  Zg)  3 


j(zp+zq) 


(  Sid  )doi 


(2.27) 


see  e.g.  Nash  and  Jordan  (1978). 

2.4  Relationship  Between  Degree-Variances  and  the  Power  Spectrum 

A  unique  and  simple  relationship  exists  between  the  spherical  harmonic  degree- 
variances  and  the  flat-earth  power  spectrum.  This  relationship,  which  says  that 
the  degree- variances  and  the  power  spectrum  are  more  or  less  the  same,  have  been 
given  e.g.  by  Dorman  and  Lewis  (1970),  in  a  form  slightly  different  than  the  one 
given  below. 

Consider  a  local  covariance  function  K(s),  where  s  is  the  distance.  The 
power  spectrum  is 


^^(u)  3  2ir7(s)  3  2tt  f"  sK(s)J.(wS)ds 
II  n  0 


while  the  corresponding  degree- variances  are  by  (2.4) 


(2.28) 


o.  -  -=M-  /"  K(ipR)  P  (cos  *)sin<frd<f> 

*  C  n  t 


(2.29) 


•I-.-vV-S 


•  *v  c* .*  v s  v  s'  s“  « 


f.f. 
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where  the  covariance  function  is  now  viewed  as  a  function  of  s  rather  than  ^  . 
However,  for  a  local  covariance  function,  obtained  e.g.  when  a  reference  field 
has  been  subtracted,  the  function  K(^R)  will  be  virtually  zero  for  large  distances. 
We  may  therefore  approximate  sinij,=  41,  and  obtain 


f  K(*R)  P„(cos  *)4>d<j,  (2.30) 

The  function  K  has  here  formally  been  extended  to  infinity.  Now  (2.30)  is  trans¬ 
formed  using  the  asymptotic  expansion  of  Legendre  functions  in  terms  of  Bessel 
functions  (Gradstein  and  Ryshik,  1965,  8.722) 


Pt(cos  *)  =  J0((2£+l)  sin  |)  ♦  0(d;2 )  =  J0(-^  4-)  (2.31) 


This  approximation  is  valid  with  high  accuracy,  even  for  relatively  low  degrees 
and  large  distances.  For  example  £=10,  ^=8. 1°  we  have 


(cos  4) )  =  0.5201 

J0(-^4O  =0.5196 


giving  only  a  relative  error  of  1%. 

By  insertion  of  (2.31)  into  (2.30)  then 


o 


£ 


/"  K(4<R)Jo{-^-  4')4»d4< 

/;  K(s)J0(^is)sds 

^  F  1ST  ^TT 


(2.32) 

(2.33) 

(2.34) 


Thus,  the  degree- variances  are  obtained  essentially  simply  by  multiplying  the 

o 

power  spectrum  by  £,  at  the  "natural"  wave  numbers 


.... 


Therefore,  results  derived  in  the  sequel  for  the  power  spectra  of  various 
density  distributions  and  masses  may  be  directly  transformed  into  “spherical" 
results,  especially  since  the  emphasis  here  is  on  local  phenomena.  As  seen  from 
the  numerical  example,  the  derived  degree- variances  may  be  sufficiently  accurate 
even  as  far  down  as  *=10. 


3.  The  Relationship  Between  Covariance  Functions  and  Mass  Distributions 

The  anomalous  gravity  field  is  generated  by  variations  in  the  density  dis¬ 
tribution  inside  the  earth.  Formally,  we  may  write  T  as  an  integral  over  the 
earth 

T  *  G  J  AP-  dn  (3.1) 

a 

where  r  is  the  distance,  G  the  gravitation  constant  and  Ap  the  density  anomaly,  i.e. 

Ap  *  p  -  pQ  (3.2) 

where  p  is  the  actual  density  and  pQ  a  normal  density  distribution,  generating 
the  normal  potential.  Opposed  to  p,  Ap  may  naturally  attain  both  positive  and 
negative  values.  When  a  reference  field  is  utilized,  the  density  anomalies  may 
be  viewed  as  being  relative  to  a  reference  density  distribution,  generating  the 
reference  field.  The  reference  distribution  need  usually  not  be  given  explicitly, 
it  suffices  simply  to  view  it  as  the  expected  "normal"  structure  of  the  earth's 
interior.  For  more  details,  see  e.g.  Forsberg  (1984). 

In  the  sequel,  the  discussion  of  density  distributions  etc.  should  thus  be 
understood  in  terms  of  density  anomalies  rather  than  just  physical  densities. 

3.1  Deterministic  Power  Spectra  of  Elementary  Masses 

Power  spectra  of  simple  anomalous  density  bodies  have  widespread  use  in  geo¬ 
physical  interpretation  and  also  play  a  role  In  understanding  the  different  types 
of  covariance  functions  frequently  used  In  geodesy.  Consider  first  a  point  mass, 
or  -  equivalently  -  a  sphere  of  constant  density.  The  gravity  effect  in  P  for  a 
sphere  of  mass  m  will  be  (cf.  Figure  2) 


vertical  DIPOLE 


VERTICAL  MASS  LINE 


'J*  "Jl  "  J*  "J* 
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g  =  Gm: 


(3.1) 


and  thus  by  (2.23) 

g(oj)  =  2w6me'uD  (3.2) 

For  the  gravity  and  potential  power  spectra  thus 


*gg(“)  -  li(“)l2  =  (2irGin)2e*2“D  -  e'2u)D  (3.3) 


$TT(ui)  *  — j-  $  (oi)  — e 

tt  '  oi2  gg  /  w2 


(3.4) 


where  is  used  also  for  "proportional  to". 

For  the  vertical  dipole  of  moment  u  we  have 
by  well-known  potential  theory 


T  =  Gp 


(3.5) 


Figure  2 

and  again  for  the  spectrum 


-u>D 


T(oi)  =  2irGue" 

In  this  case  the  power  spectra  have  more  power  at  higher  wave  numbers: 


(3.6) 


2  -2aiD 

gg ' 

-2ujD 

TT  ' 

e 

(3.7) 


■-1 


v.K 


r-’-l 


•’>0 


‘  (3.8) 
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The  vertical  mass  line  has  less  power  at  higher  wave  numbers.  We  have 

g  s  Gk  r  Adz  =■  Gk~—  (3.9) 

0  r  rtop 

where  <  is  the  line  mass  density.  Thus 


99  ’ 

1  -2u)D 

2  ® 

(3.10) 

1  -2u>D 

(3.11) 

TT  * 

4  ® 

point  mass 


Conventionally  in  potential  field  analysis 
(gravity  and  magnetics),  power  spectra  are 
cos  \  plotted  (single)  logarithmic.  In  this  case 

( da )  v  \mosi  hne  the  point  mass  power  spectrum  will  be  a  straight 

N.  '  ,  line,  the  slope  of  which  is  a  direct  measure 

^\V  of  the  depth  D  (Figure  3).  This  forms  the 

S*  base  spectral  methods  in  geophysical  inver- 

jdipoie  sion,  especially  used  in  magnetics.  (In  the 

- —  <jj  magnetic  case,  gravity  results  may  be  transfer¬ 
red  by  using  Poissons  formula  for  the  magnetic 
Figure  3  potential:  U*-m  *  ?T,  where  m  is  the  magnetization 

vector  of  the  body  in  question.  This  implies 
e.g.  that  the  gravity  point  mass  spectrum  corresponds  to  a  magnetic  vertical  dipole 
line). 

It  is  of  interest  also  to  consider  the  change  in  shape  of  the  power  spectra 
when  the  dimensions  of  the  "anomalous  body"  are  not  negligible.  To  include  "width" 

and  "thickness"  of  the  body,  consider  a  cylinder 


Figure  3 


VERTICAL  CYLINDER 
s  P 


of  constant  density  p,  diameter  2a  and  thickness 
T  (Figure  4).  The  gravity  effect  is  one  of 
the  less  trivial  geophysical  elementary  bodies 
(the  "volcanic  plug").  Using  (3.9)  the  attrac¬ 
tion  at  P  may  be  written 


Figure  4 


g(s)  =  Gp  /  (-±-  -  -^-)2ira'da' 
o  rt  rb 


(3.12) 


r\  -  (s  -  a1)2  +  D2, 


r2  =  (s  -  a1)2  +  (D  +  T)2 


To  get  the  spectrum,  (3.12)  is  inserted  into  the  Hankel  transform  integral  (2.19): 


gU)  =  2ng(to)  =  4it2Go  J  /  (~-  -  -7-)  a'da'J  (u>s)sds 

0  0  rt  rb  0 


(3.13) 


which  by  interchange  of  integration  order  and  some  evaluations  results  in  (Petersen, 
1978) 


g(»)  ■  2,0.  •  e-“°  •  iC  •  iWial, 
3  011  ^00 


m  =  ira2T 


(3.14) 


This  is  seen  to  be  the  point  mass  spectrum  (3.2)  multiplied  by  a  thickness  and 
a  width  factor  (Jt  is  the  Bessel  function  of  the  order  1).  These  factors  are 
shown  in  Figure  5  for  unit  radius  and  thickness.  The  factorization  is  typical 

for  the  spectra  of  "simple"  bodies.  A 
similar  formula  is  obtained  e.g.  for  the 

1  ^  rectangular  prism,  see  e.g.  Forsberg  (1984). 

\\widlh 

\\  Due  to  the  factorization,  the  power  spec¬ 
s'  trum  becomes  equally  simple: 


Jhickness 


8  10 


Figure  5  When  T  «  D  and  a  «  D,  the  factors  may  be 

neglected,  or,  in  other  words,  the  point  mass  approximation  is  adequate. 

3.2  Statistical  Mass  Distributions 

The  spectra  of  the  "simple"  bodies  of  the  last  section  are  encountered  also 
when  considering  independent  statistical  ensembles  of  the  same  bodies.  This  is 
due  to  a  fundamental  theorem  of  statistical  mechanics,  stating  that  the  mathematical 


■L- ^  (vS'! 


vv. 
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expectation  of  the  "total"  power  spectrum  equals  the  ensemble  average  of  the  "Indi¬ 
vidual"  power  spectra. 

Assume 


4  M 


$P(u>.  Pj , 


Pk} 


(3.16) 


Is  the  power  spectrum  of  an  individual  body,  characterized  by  parameters  p  ,  ...»  pk, 
describing  e.g.  depth,  width,  density  etc.  (radial  symmetry  has,  for  simplicity, 
been  assumed).  The  distribution  of  the  bodies  is  described  by  an  ensemble  joint 
frequency  distribution  f (pi ,  ....  pk).  Then  the  total  power  spectrum  $  will  be 
the  ensemble  average  of  : 


P  P 

4>U)  *  E{$  («)}  ■  J  •••  /  <p  (w,  . pk)  f(Pj,  ...  ,  pk)  dpj  ...  dpk 

Pj  Pk  (3.17) 

(Spector  and  Grant,  1970). 

If  the  parameters  are  assumed  to  vary  independently,  then  the  joint  frequency 
distribution  will  be  a  product  f(plf  ...»  pk)  *  fl(P1)  fk(Pk)*  0f  special 
interest  is  the  case  where  the  parameters  are  assumed  to  have  a  white  noise  dis- 

p 

tribution.  In  this  case  $(u>)  ■  <f>  (w),  and  we  have  immediately  the  following  spectra 
for  a  “random"  layer  at  depth  D: 


White  noise  density  layer: 

.  I  \  -2(uD 

♦gg(“)  -  e 

(3.18) 

White  noise  dipole  layer: 

/  »  2~”2u>D 

*ggW)  ' 

(3.19) 

White  noise  "mass  lines": 

/  \  1  -2u)D 

^ggU)  *  H7  e 

(3.20) 

The  white  noise  density  layer  may  be  interpreted  physically  as  random  undulations 
in  the  depth  of  the  interface  between  two  layers  of  constant  density,  e.g.  a  sedi¬ 
ment/basement  interface. 


0 


Figure  6 


For  an  example  with  two  parameters,  consider  an 
earth  model  of  point  masses  of  white-noise  random 
mass  and  random  depth  below  a  layer  at  D  (Figure  6) 
For  this  "all  white"  assemblage  of  mass  points 

.  /  \  r  -2u)D  jni  1  -2wD  i  i  r 

OggU)  '  J  6  °D  '  W  6  (3.< 


and  similarly  for  an  "all-white"  dipole  distribution 


♦gg<“)  '  “e 


-2u>D 


(3.22) 


This  may  be  generalized  to  the  so-called  stationary  thin  layered  earth  model, 
where  the  density  anomaly  distribution  is  assumed  to  be  described  by  a  covariance 
function  and  thus  power  spectrum  $  .  If  the  earth  is  assumed  to  consist  of  a 

number  of  thin  layers  (e.g.  a  sedimentary  sequence),  then  each  layer  will  yield 
a  partial  power  spectrum  of  form 


<>» _ (oo)  =  (2ttGaz)2  0  (a),  d)e"2wd  (3.23) 

gg  pp 

where  d  is  the  depth  and  az  the  thickness. 

If  the  layers  are  similar  and  independent,  the  total  power  spectrum  thus 
will  be 


V»>  -  ♦„„<•> ;  e'2"°  <3-24> 

yielding  an  example  of  a  one-to-one  relationship  between  gravity  and  density  co- 
variance  functions. 

The  above  examples  all  have  in  common  the  factor  e“2cuD,  containing  the  depth 
information.  This  factor  is  the  prime  "objective”  of  the  widespread  statistical 
inversion  techniques  of  geophysical  inversion  applied  successfully  for  more  than 
a  decade,  especially  for  aeromagnetic  data.  In  the  method,  depths  to  the  top 


and  base  of  magnetized  layers  may  be  found  (at  times  also  width  parameters  are 
solved  for),  using  the  shape  of  the  power  spectrum,  computed  for  a  given  area. 
Layers  at  various  depth's  typically  show  up  as  more  or  less  straight  line  segments 
in  the  power  spectrum,  when  plotted  logarithmic,  and  depths  are  derived  from  the 
slope  values.  Examples  may  be  found  e.g.  in  Spector  and  Grant  (1970). 

An  impressive  example  of  a  two-layer  magnetic  problem  may  be  seen  in  the 
global  degree  variances  of  the  magnetic  field  (Figure  7).  The  field  from  the 
earth's  core  is  seen  to  dominate  up  to  z=13,  after  which  the  "shallow"  crustal 
field  takes  over.  The  core  radius  Inferred  from  the  slope  of  the  low  degree 
variances,  yields  a  value  only  a  few  hundred  kilometers  deeper  than  the  seis- 
mically  determined  value  (3485  km). 
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Figure  7  Degree-variances 
for  the  earth's  magnetic  field. 
MAGSAT  solution,  complete 
to  degree  and  order  23  (after 
Langel ) 


While  the  statistical  approaches  have  proved  very  successful  for  magnetic 
data,  they  have  been  less  applicable  to  gravity  interpretation.  This  is  partly 
due  to  the  (usually)  less  dense  gravity  coverage,  and,  especially,  to  the  dif¬ 
ferent  nature  of  the  sources:  density  anomaly  variations  tend  to  be  more  "spread 
out"  and  less  random  than  changes  in  magnetization  and  susceptibility.  The 


typical  “clear-cut"  geophysical  problem  of  unmagnetic  sediments  overlying  a  highly 
magnetized  methamorphic  basement  has  e.g.  usually  a  much  less  clear  gravity  paral¬ 
lel,  intra-sedimentary  density  anomalies  being  most  often  non-negligible. 

However,  the  value  of  the  statistical  interpretation  for  gravity  is  yet  high: 
it  provides  an  efficient  means  for  a  deeper  understanding  of  covariance  functions, 
as  outlined  in  the  sequel. 

3.3.  Relationship  to  Traditional  Covariance  Functions 

Most  frequently  simple,  analytical  covariance  functions  have  been  used  for 
e.g.  gravity  field  modelling  by  collocation.  These  functions  are  characterized 
by  a  few  "free"  parameters,  which  may  be  adjusted  to  fit  the  "essential"  parameters 
of  a  local  covariance  function,  such  as  (gravity)  variance  C0,  correlation  length 
x :.  and  horizontal  gravity  gradient  variance  G0,  for  a  discussion  see  Moritz  (1980, 
ch.  22). 

For  the  planar  case  many  different  types  of  analytical  covariance  functions 
have  been  suggested:  e.g.  exponential,  Gaussian  and  especially  models  of  form 


C(s) 


C0 

(l+(s/D*)2)m 


(3.25) 


where  O'  and  m  are  constants.  However,  of  these  simple  models,  only  two  allow 
the  derived  spatial  covariances  C(s,  z^  z2)  to  have  correspondingly  simple  ex¬ 
pressions:  namely  (3.25)  for  m=l/2  and  m=3/2,  see  again  Moritz  (1980).  The  co- 
variance  functions 


C(s) 


Cp 

(H^-)2)*5 


and 


C(s) 


Co 

(H^r)2)^ 


(3.26) 


(3.27) 


have  been  termed  the  reciprocal  distance  and  the  Poisson  covariance  functions 
respectively.  The  associated  power  spectra  follow  from  (2.22)  and  (2.23): 
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1  —  D 1 

Reciprocal  distance:  *  2TrC0D'^-e*w  (3.28) 

Poisson:  *  2irC0D'2e'“D  (3.29) 

Thus,  by  (3.22)  the  reciprocal  distance  covariance  function  may  be  interpreted 
as  generated  by  an  "all-white"  assemblage  of  mass  points  from  depth  D=y)‘,  or 
in  other  words  -  a  stationary,  thin-layered  earth  with  a  white  noise  density  anomaly 
distribution  below  depth  D.  The  Poisson  covariance  function  may  similarly  be  in¬ 
terpreted  as  associated  with  a  white  noise  density  layer  at  depth  D=y)'. 

For  spherical  covariance  models  a  number  of  expressions  have  been  used  or 
suggested.  Generally,  the  covariance  models  have  been  obtained  from  selected 
models  of  the  degree- variances,  where  the  spatial  covariances  may  be  evaluated 
by  analytical  expressions. 

At  sea  level  the  general  expression  for  the  potential  covariance  function 
may  be  written 


K(*)  *  l  at  $)*+1  Pfc(cos  *) 
1  2 


(3.30) 


where  R^  is  the  radius  of  the  Bjerhammar  sphere  and  R > an  earth  radius.  The 
depth  to  the  Bjerhammar  sphere  is  directly  related  to  the  depth  parameter  D  of 
the  planar  covariance  functions:  When  R-R^  «  R  then 

(^)*’+1  »  exp(2(i+l)  log  A))  =  e'2u)D,  u  «  D  =  R-R.  (3.31) 

K  R  u 


Thus  for  large  i,  the  asymptotic  degree  variances  for  the  reciprocal  distance 
and  Poisson  covariance  functions  must  be  of  the  form 


Reciprocal  distance: 


a  „  Jj.  0*+1  (o  =  $) 


(3.32) 


1  l+i 
°t  *  7  0 


Poisson: 


(3.33) 


This  follows  directly  from  the  relationship  (2.34)  between  the  power  spectrum 
and  the  degree  variances.  Closed  formulas  for  the  spherical  reciprocal  distance 
and  Poisson  covariance  functions  are  found  in  Moritz  (1980,  ch.  23). 

An  Important  additional  class  of  degree-variance  models  is  the 


Logarithmic:  -  jy  ai+l 


(3.34) 


The  earth's  gravity  field  is  known  to  adhere  to  such  a  power  law  on  a  global  scale 
Well  known  empirical  logarithmic  models  are  "Kaula's  Rule" 


.  n  7  in-io  2i+l 
0„  -  0.7  *  10  -r-TT~ 


(3.35) 


(see  e.g.  Phillips  and  Lambeck,  1980)  and  the  "Tscherning-Rapo"  model 


=  4.4  •  10 


-10 


1 


U-l)U-2)(*+24) 


.99962 


z+i 


(3.36) 


(Tscherning  and  Rapp,  1974),  which  corresponds  to  a  Bjerhammar  sphere  depth  of 
1.2  km.  From  (3.34)  follows  that  the  asymptotic  form  of  the  logarithmic  power 
spectrum  is 


<P  ((a)  )  *• 

vggv  ' 


Z1 


■2(1)0 


(3.37) 


which  is  seen  to  be  the  "white  mass  lines"  power  spectrum.  A  more  reasonable 
physical  model  is  obtained  by  the  stationary,  thin-earth  assumption:  then  from 
(3.24)  the  density  spectrum  must  be  of  the  form 


J. 

01 


(3.38) 


which  Is  a  plausible  "red  noise"  spectrum  with  high  power  in  the  low  wave  numbers. 
However,  such  a  noise  model  globally  implies  stresses  in  the  Interior  of  the  earth 
that  are  bigger  than  the  stress  which  e.g.  the  mantle  is  currently  thought  to 


be  able  to  sustain:  therefore  "multilayer"  models  are  more  adequate  for  obtaining 
realistic  statistical  density  models.  Naturally  such  models  in  addition  may  provide 
a  better  fit  to  empirical  data,  since  more  free  parameters  will  be  available  to 
"tune"  the  covariance  model. 

Heller's  attenuated  white  noise  (AWN)  model  is  an  example  of  this  kind. 

In  the  AWN  model  the  global  degree-variances  are  modelled  using  5  density  shells 
at  various  depths,  where  the  potential  originating  from  each  shell  is  assumed 
to  be  white  noise  at  the  shelf  itself,  i.e. 

a  ”  l  ai(2l+l)(-^^-)"+1  (3.39) 

1  i=l  1  R 


where  a.  and  D  are  constants  (Jordan,  1978).  For  each  shell,  the  associated 

I  i) 

planar  power  spectrum  will  be  of  form 
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w2e 


-2wDi 


(3.40) 


and  thus  the  AWN-model  may  be  interpreted  as  a  white  noise  dipole  layer  model . 

The  various  simple  covariance  models  are  listed  in  Table  1,  showing  the  un 
normalized  forms  of  the  gravity  power  spectrum  $  (w),  covariance  function  C(s) 
and  degree-variances  using  (2.34) 


Table  1 


Density  Model 

*qq(">) 

C(s) 

Asymtotic 

Covariance  type 
(Moritz) 

white  mass  lines 

1  -2d)D 

2  e 

-log(D'-t-r) 

1 

u 

"logarithmic" 

"all-white"  masses 

1  -2wD 
—  e 

0) 

I 

r 

1 

u 

"reciprocal 

distance" 

white  mass  layer 

-2wD 

e 

D' 

1 

l 

"Poisson" 

"all-white"  dipoles 

-2wD 

we 

3D'2  ..  1 
r^  7T 

const. 

white  dipole  layer 

2  -2wD 
w^e 

15D ' 3  3D' 

-fr~  -  yr 

i 

"AWN" 

Note:  O'  *  2D,  r  =V  D'^+s2 


The  covariance  function  expressions  have  been  given  by  (3.26)  and  (3.27) 
for  the  "reciprocal  distance"  and  "Poisson",  respectively.  The  remaining  co- 
variance  functions  have  been  derived  utilizing  the  relationships 

w  •  •« — ►  differentiation  ^r 
$gg/u>  «« — ►  integration  after  D' 

between  frequency  domain  and  space  domain.  Thus,  C(s)  one  step  downwards  in  the 
table  is  obtained  by  differentiation,  C(s)  one  step  upwards  is  obtained  by  inte¬ 
gration. 

Correspondingly,  since  for  the  potential  power  spectrum 


the  form  of  the  potential  covariance  function  K(s)  may  be  found  going  two  steps 
up  in  the  table,  and  similarly  second-order  gradient  covariances  may  be  found 
going  two  steps  down. 

When  going  "upwards"  in  the  table,  the  covariance  function  C  is  seen  to  be 
less  and  less  "sharp",  and  from  the  "logarithmic"  level  (and  upwards)  it  will  be 
singular,  going  to  minus  infinity  for  large  r.  This  singularity  is  a  consequence 
of  the  inadequacy  of  the  planar  approximation  for  low  wave  numbers  .  The  large 
power  at  low  wave  numbers  (^r  in  the  "logarithmic"  case)  will  not  be  found  on 
actual  power  spectra,  since  the  use  of  a  spherical  harmonic  reference  field  will 
remove  the  power  at  the  lowest  wave  numbers.  For  all  the  simple  power  spectra 
of  Table  1,  such  suppression  of  low  wave  numbers  will  result  in  well-behaved  co- 
variance  functions  C.  In  principle,  the  necessary  suppression  of  low  wave  numbers 
is  corresponding  to  the  omission  of  the  lowest  degree-variances  (e.g.  z-0  and  1) 
in  the  spherical  covariance  expression. 

The  suppression  of  the  low  wave  numbers  will  be  treated  in  more  detail  later 
in  this  report.  First,  however,  attention  will  be  given  to  the  covariance  expres¬ 
sions  C(s).  Their  simple  forms  allow  a  straightforward  interpretation  of  collocation 
in  terms  of  "generalized  point  mass"  modelling  under  certain  conditions. 


4.  Collocation  as  Generalized  Point  Hass  Hodelllng 

4.1  Covariances  as  Gravimetric  Effects 

The  covariance  expression  C(s)  of  Table  1  may  be  viewed  as  gravity  effects 

for  points  at  the  reference  level.  An  example:  consider  two  points  P  and  Q  with 

distance  S.  Then  the  gravity  covariance  C(P,  Q) 

for  a  white  noise  density  layer  at  depth  D 

poisson  (the  Poisson  covariance)  corresponds  to  the 

p  n 

u .... 1 1 .. T. ■■*■■■/.  ■  r  gravity  effect  in  Q  from  a  point  mass  at  depth 

2D  below  P: 


rec.dtsl. 


P  Q 


logon  Imic 


Figure  8 


“Poisson"  C(P,  Q)  -  agMQ) 
cf.  Figure  8. 

Similarly  the  "reciprocal  distance"  covariance 
function  ("all  white"  masses)  corresponds  to  the 
mass  line  effect,  while  the  logarithmic  covariance 
may  be  viewed  as  generated  by  a  rather  strange 
mass  body:  a  "wedge",  i.e.  a  mass  line  distribu¬ 
tion  with  line  mass  density  <  increasing  linearly 
downwards  (Figure  8). 

Thus,  the  "successive"  covariance  functions 
of  Table  1  are  generated  by  the  following  "general¬ 
ized  point  masses": 


"logarithmic" 
"rec.  dist." 
"Poisson" 


mass  wedge 
mass  line 
point  mass 
dipole 
quadropole 


(these  "deterministic"  mass  distributions  should, 
of  course,  not  be  mixed  with  the  statistical 
density  models  of  the  previous  section). 
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In  fact,  an^  covariance  function  C(s)  may  be  interpreted  as  some  line  mass 
effect: 

Assume  C  to  be  the  gravity  effect  of  a  mass  line  below  P  of  line  density 
x(t),  where  t*-z  is  the  depth  coordinate  (Figure  9).  Now, 

C(s)  =  G  /  tdt  (r  =  </sz+tz)  (4.1) 

o 

and  by  taking  the  Fourier  transform  and  interchanging 
the  order  of  integration 

4>  (a)  =  ZrrG  J  *(t)  e"utdt  (4.2) 

UU  0 

Again  (2.3)  has  been  used.  Thus,  the  power  spectrum 
4>gg  is  given  as  the  Laplace  transform  of  k.  Conversely, 
given  C,  «  may  be  found  by  an  inverse  Laplace  transform 
of  the  associated  power  spectrum.  Due  to  the  properties 
of  the  Laplace  transform,  the  transition  C« — is  in 
principle  unique  and  exists  for  all  reasonable  covariance  functions. 

The  gravimetric  interpretation  also  holds  for  covariances  aloft  under  certain 
conditions.  If  P  is  at  reference  level,  but  Q  at  altitude  z,  then  by  (2.27): 


/  \  Cl*) 

✓  \ 

/  V 

P 


wp 


Figure  9 


C(P,  Q)  =  ^  e  aiZJ(J(sw)du 

o  ** 


By  inserting  (4.2)  and  interchanging  the  order  of  integration  then 


C(P,  Q)  *  Gf  *(t)  J  we”w(t+z)Jn(su.)da,dt 


=  G  J"  - ^ — sr  (t+z)dt 

o  [sz+(t+z)z]  ^ 


(4.3) 


(4.4) 

(4.5) 


which  is  simply  the  gravity  effect  at  Q.  If  both  P  and  Q  are  at  altitude,  then 
z  in  (4.5)  should  be  replaced  by  zp+Zp.  The  "gravity  interpretation"  of  a 
covariance  function  C(s)  thus  allows  an  easy  extension  to  the  spatial  covariance 
C(s,  Zp,  Zg). 
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4.2  Application  to  Least  Squares  Collocation 

In  the  statistical  formulation  of  least  squares  collocation  a  signal  s  is 
estimated  from  a  vector  of  measurements  x,  the  estimate  being  given  by  the  well- 
known  formula 


CsxCxxx 


(4.6) 


The  cross-  and  auto-covariances  C$x  and  Cxx  are  obtained  by  covariance  propogation 

!c»>r  w'--  •>  {4-7) 

<Cxx’ij  ■  LP1LPjK<’ •  ■’  (4'8) 


where  the  signal  s=Ls(T)  and  observations  x^=Lpi. (T)  are  viewed  as  linear  functionals 
(the  dots  in  (4.7)  and  (4.8)  indicate  that  the  functionals  should  be  applied  on 
separate  variables  of  the  potential  covariance  function  K(P,  Q)). 

For  practical  solutions,  s  is  obtained  as  a  linear  combination 


5  *  |  w, 


(4.9) 


with  coefficients  ai  given  as  the  solution  to  the  "normal  equations" 


{{CXX>U  a1  ’  xj  (4-10) 

Consider  now  the  simple  case  of  the  observations  being  gravity  anomalies 
at  a  set  of  points  P^ ,  and  assume  the  Poisson  covariance  function  to  be  used. 
Then,  with  the  interpretation  of  covariances  in  terms  of  gravity  effects,  the 
covariances  may,  as  earlier  mentioned,  be  written  as 


(Cxx>,j  ■ 


(4.11) 
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where  qg1  represents  the  gravity  effect  of  a  mass  point  of  unit  mass,  located 
below  P^  at  a  depth  D‘  determined  by  the  chosen  correlation  length  of  the  Poisson 
covariance  function  used,  i.e.  D'  *  1.30  x^,  as  may  easily  be  verified. 

Then  (4.10)  may  be  written  as 

l  agVj)  ai  =  x j  (4.12) 

which  is  nothing  but  the  inversion  equation  of  point  mass  modelling,  with  point 

masses  located  below  all  points  P..  Thus  collocation  and  point  mass  modelling 
•  J 

are  identical  when  the  Poisson  covariance  function  is  used,  and  signal  estimates 
(4.9)  are  simply  obtained  as  the  relevant  computed  gravity  field  quantity  from 
the  point  masses,  the  mass  point  below  the  i'th  gravity  point  having  mass  a^ . 

For  other  types  of  covariance  functions,  collocation  may  similarly  be  viewed 
as  generalized  point  mass  modelling,  where  point  masses  are  then  replaced  by  the 
other  "elementary"  mass  bodies  (cf.  Section  4.1).  In  e.g.  the  AWN  case,  colloca¬ 
tion  thus  corresponds  to  solving  for  a  number  of  quadropole  moments  while  e.g. 
the  logarithmic  covariance  functions  correspond  to  “mass  wedge"  modelling. 

The  above  only  hold  when  exclusively  gravity  data  are  used  as  observations. 

For  heterogenous  data,  e.g.  a  mixture  of  gravity  and  geoid  observations,  the  simple 
scheme  must  be  modified  so  that  the  type  of  elementary  mass  model  to  be  used  below 
an  observation  point  will  be  dependent  on  the  type  of  gravity  field  data  given 
at  the  point. 

Consider  again  the  Poisson  covariance  function,  and  assume  a  mixture  of  gravity 
and  geoid  data  to  be  given  at  the  reference  plane,  e.g.  representing  available 
gravity  field  data  in  an  oceanic  area.  Then  the  covariances  in  Cxx  between  gravity 
stations  are  identical  to  point  mass  effects.  Between  points  with  geoid  observa¬ 
tions,  the  covariances  will  be  (essentially)  K(<p ) ,  which  will  be  of  form  -log(D'+r), 
cf.  Table  1,  which  is  the  form  of  the  geoid  effect  of  the  vertical  mass  line  with 
top  at  depth  D‘.  Similarly,  geoid-gravity  cross  covariances  may  be  interpreted 
as  the  gravity  effect  of  the  vertical  mass  line,  or  -  equivalently  -  the  geoid 
effect  of  the  point  mass. 

In  this  case  thus  (4.10)  may  again  be  viewed  as  an  inversion  equation,  where 
the  unknowns  to  be  solved  for  will  be  for 

Gravity  stations:  point  mass  values 


Geoid  points  :  line  density  values 
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as  shown  schematically  in  Figure  10.  Other  gravity  field  data  types  have  similar 
associated  "Inversion  masses":  second-order  gradients,  e.g.  Tzz* corresponds  to 
a  dipole  point,  while  deflections  of  the  vertical  like  gravity  corresponds  to 
the  point  mass. 


pornon  X  GRAVITY 

m  OEOID 


•  •  I  •  I  •  D' 


Figure  10  Collocation  as  generalized  point  mass  modelling.  With  the  Poisson 
covariance  function  ("white  noise  density  at  depth  D"),  collocation 
corresponds  to  solving  for  point  masses  and  line  mass  densities  at 
depth  D' =20,  satisfying  the  given  observation  data. 

The  view  of  collocation  as  a  generalized  point  mass  modelling  technique  is 
of  course  primarily  of  interest  as  a  tool  for  a  better  clarification  of  what  is 
"going  on"  when  using  least  squares  collocation.  Naturally,  the  practical  work 
remains  the  same,  namely  to  solve  a  system  of  linear  equations  containing  one 
equation  for  each  observation  point  and  predict  "signals"  by  summing  up  the  solu¬ 
tion  coefficients  times  the  relevant  cross-covariances. 

While  collocation  in  general  may  be  viewed  as  generalized  point  mass  model¬ 
ing,  traditional  mass  modelling  approaches  are  not  necessarily  corresponding  to 
collocation.  Only  when  the  number  of  unknowns  equals  the  number  of  observations 
(and  each  "mass  point"  is  situated  Immediately  below  an  observation  point)  will 
the  result  of  point  mass  modelling  be  a  collocation  solution:  point  mass  modelling 
of  gravity  data  corresponds  to  the  use  of  Poisson's  covariance  function,  dipole 


modelling  to  the  "all-white  dipoles"  covariance  function  (one  step  down  in  Table  1) 
etc.  Since  the  earth's  actual  gravity  field  variations  tend  to  be  best  described 
by  logarithmic  covariance  functions,  "optimum"  results  should  thus  be  obtained 
with  some  "mass  wedge"  type  inversion. 

5.  Application  of  Simple  Covariance  Models  for  Local  Covariance  Functions 

After  the  "side  jump"  of  the  last  section  to  collocation,  emphasis  will  now 
return  to  the  statistical  density  distributions  and  the  associated  covariance 
functions. 

5.1  Use  of  Multi -Layer  Models 

The  simple  covariance  functions  of  the  previous  sections  may  be  used  to  build 
good  approximations  to  empirical  covariance  functions  through  the  use  of  multi¬ 
layer  models. 

Consider  the  simple  case  of  a  number  of  white  noise  density  layers,  representing 
e.g.  interfaces  at  a  layered  earth  model.  If  the  undulations  of  the  interfaces 
are  assumed  statistically  independent  (Tike  the  "thin-layered  stationary  earth 
models"),  then  the  total  covariance  function  will  be  the  sum  of  the  Individual 
(Poisson)  covariance  functions 

C(s)  -  I  Us)  -Ig-frl.  r.  -v'F+TCTF-  (5.1) 

1  1  i  Em1  ri  1  1 

as  illustrated  in  Figure  11  (m^  is  the  variance  of  each  mass  plane  density  dis¬ 
tribution). 


Figure  11  Two-layer  stationary  density  model  with  associated  power  spectrum  and 
covariance  function. 
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6y  adjusting  the  "free"  parameters  (D. ,  m^)  good  approximations  to  empirical  co- 
variance  f uncations  may  be  obtained.*) 

The  depth  and  density  parameters  may  be  chosen  based  on  geological  considera¬ 
tions  or  -  usually  more  realistic  -  based  on  the  shape  of  the  computed  empirical 
power  spectra.  If  the  gravity  field  sources  at  a  certain  depth  range  tend  to 
behave  like  a  white  noise  layer,  then  the  power  spectrum  tend  to  show  straight- 
line  segments  at  a  wavelength  band  corresponding  to  the  depth  range,  when  plotted 
logarithmic. 

An  actual  data  example  (Ohio  area)  is  shown  in  Figure  12.  In  this  example 
a  reasonably  good  covariance  fit  is  obtained  using  a  two  layer  model,  with  depths 
determined  from  the  shape  of  the  two  fitted  lines.  Below  0.5  cycle/degree,  cor¬ 
responding  to  spherical  harmonic  degree  180,  the  model  should  not  fit,  since  har¬ 
monics  below  180  have  been  attempted  to  be  removed  by  utilization  of  a  180x180 
reference  f igld  (  the  apparent  existence  of  power  below  0.5  cycle/degree  is  due 
to  errors  in  this  reference  field). 

In  the  Ohio  example,  the  two  depth  values  represent  the  average  influence 
of  shallow  and  deep  sources.  Considering  the  geology  (primarily  a  paleozoic  sedi¬ 
mentary  basin  of  thickness  of  order-of -magnitude  DJ,  the  two  covariance  model 
constituents  may  be  ascribed  as  originating  primarily  from  sedimentary  sources 
and  deep  crustal  sources,  (e.g.  undulations  of  the  "conrad"  discontinuity)  respec¬ 
tively.  Such  an  interpretation  should,  however,  be  taken  with  great  caution. 

The  area  is  so  big  that  there  is  no  reason  to  expect  the  geology  to  be  so  "constant" 
that  the  layering  of  the  crust  "shines  through"  to  the  estimated  power  spectrum, 
and  also  the  shallow  depth  parameter  is  affected  by  important  error  sources  such 
as  data  noise  and  aliasing. 

With  the  two-component  description  of  the  power  spectrum,  it  is  easy,  by 
simple  graphical  techniques,  to  subdivide  the  gravity  field  variation  into  two 
parts,  representing  each  source. 

Each  component  is  a  Poisson  covariance  function  of  the  form 


C(s) 


r  u 

0  [s*+0,ajfc 


(0**20) 


(5.2) 


7) - 

'Note  that  with  e.g.  simple  Poisson  covariance  functions,  C(s)  can  never  be  negative. 
However,  for  practical  applications  the  covariance  functions  need  to  be  modified 
by  removing  power  at  the  lowest  wave  numbers,  to  avoid  the  singularities  in  the 
elementary  functions  caused  by  the  deficiency  of  the  planar  approximation  for  low 
wave  numbers.  After  this  "removal"  negative  covariances  may  be  obtained  as  well. 
More  details  in  the  next  section. 
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Figure  12  Empirical  power  spectrum  for  residual  free-air  anomalies  with  respect 
to  Rapp's  180x180  spherical  harmonic  expansion.  Area  bounded  by 
latitude  38°-42°  N  and  longitude  85°  to  81°  W.  Data  points  (shown 
with  dots)  obtained  from  4'x4‘  gridded  gravity  data  using  a  2-dimen¬ 
sional  FFT  and  radial  smoothing. 


which  as  earlier  mentioned  has  the  power  spectrum 


»„(•)  ■  e'"D’ 


Looking  at  Figure  12,  the  ratio  between  the  two  components  at  DC  U=0)  is  seen 
to  be  around  23dB  which  by  (5.3)  results  in  a  gravity  variance  ratio  of  0.099  between 
the  shallow  and  deep  components,  respectively.  From  the  total  gravity  variance 
(153  mgal2)  the  r.m.s.  variation  of  the  shallow  and  deep  components  then  becomes 
3.7  and  11.8  mgal  respectively. 

5.2  Modifications  to  Eliminate  the  Influence  of  Low  Freauencies 


If  one  tried  to  estimate  geoid  variances  for  the  Ohio  example,  one  would 
find  that  the  derived  variances  turned  out  to  be  infinite.*)  What  is  wrong  is 
that  the  “elementary"  Poisson  covariance  function  contains  too  much  power  at  the 
lowest  frequencies  near  DC,  where  the  planar  approximation  is  not  valid. 


"^The  associated  potential  covariance  function  K(s)  will  be  of  form  -log(D’+r), 
cf.  Table  1. 
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This  is  not  the  case  for  all  of  the  covariance  models  of  Table  1.  For  in¬ 
stance  the  AWN  model  has  well-behaved  geoid  variances  even  in  the  planar 
case.  On  the  other  hand,  however,  the  logarithmic  covariance  function  is  even 
for  gravity  anomalies  singular,  since  the  derived  covariance  function  C(s)  does 
not  go  to  zero  for  large  s,  cf.  Table  1. 

There  is  therefore  a  need  for  modification  of  the  simple  covariance  models, 
so  that  the  power  at  the  low  frequencies  is  suppressed.  This  is  necessary  not 
only  due  to  the  singularities  induced  by  the  flat-earth  approximation  -  it  is 
always  needed  for  local  covariance  functions  when  a  spherical  harmonic  reference 
field  has  been  removed,  irrespectively  of  whether  a  flat-earth  or  a  round-earth 
formulation  is  used. 

If  the  spherical  harmonic  reference  field  was  assumed  to  be  error-free,  the 
natural  modification  would  be  to  truncate  the  "elementary"  power  spectra  like 
(5.3)  below  the  wave  number  corresponding  to  the  maximal  degree  of  the  spherical 
harmonic  expansion.  The  covariance  function  C(s)  is  then  derived  by  Hankel 
transformation  of  the  truncated  power  spectrum,  e.g.  for  a  Poisson  covariance 

0  OJ  OJj. 

"  ,5-4) 

2irC  D  ze  u>  >  tii a. 

0  t 

where  the  truncation  wave  number  for  a  reference  field  complete  to  spherical 
harmonic  degree  %ax  should  be  u>t  =  . 

With  this  approach  there  are  drawbacks:  First,  no  simple  analytical  expres¬ 
sions  exist  for  the  covariance  functions,  necessitating  the  use  of  e.g.  numerical 
integration  techniques,  and  secondly  for  practical  applications  some  energy 
should  be  left  below  o^,  to  take  into  account  errors  in  the  spherical  harmonic 
reference  field  (these  errors  are  believed  to  be  of  the  same  order  of  magnitude 
as  the  degree-variances  themselves  for  the  higher  degrees  in  existing  180x180 
solutions,  see  e.g.  Rapp,  1981). 

Relatively  simple  analytical  expressions  may 
p  s  q  be  obtained  by  an  approximative  alternative 

........... approach  based  on  a  pseudo-isostatic  formula¬ 
tion.  Assume  the  density  anomalies  to  be  de- 

density  *  _ ir 

layer  o  scribed  by  a  white-noise  layer  model  at  depth 

0,  and  assume  that  this  layer  is  perfectly 
Isostatically  compensated  at  a  deeper  level 
^ -  oc  0C>  The  power  spectrum  of  this  "compensated 
Poisson  model"  will  be  of  form 
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Figure  13 
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4.33(0.)  “  ar«(e-“°  -  e"wDc)2 


(5.5) 


where  a  is  a  constant.  The  formula  follows  simply  from  (3.2),  since  the  spec¬ 
trum  due  to  the  ensemble  averaging  theorem  will  have  the  same  form  as  the  spec¬ 
trum  of  a  positive/negative  mass  point  pair.  The  spectrum  (5.5)  will  be  a  linear 
combination  of  three  simple  poisson  spectra  of  depths  D,  Dc  and  Dm=J*(D+Dc): 


4.33(0.)  =  2^a(e"2w0  +  e'2uDc  -  2e  _2“D"') 


The  power  spectrum  may  adequately  be  written 


♦  U)  =  2irae-2aiD{l  -  e“wT)2,  T=D-Dc 


to  show  more  directly  how  the  introduction  of  the  compensating  layer  corresponds 
to  suppression  of  power  at  low  wave  numbers  (Figure  14). 

The  choice  of  compensation  depth  will  be  dependent 
LOG  x  on  *max  for  the  spherical  harmonic  reference  ex- 

™  pansion.  For  e.g.  a  180x  180  field,  u)t  may 

yPo/non  arbitrarily  be  put  at  the  point  where  the  "com- 

3.dB  .  compensated  pensated"  power  has  dropped  3dB  (i.e.  half) 

pomt  f. 

/;  below  the  simple  Poisson  expression,  i.e. 


Figure  14 


-2wtD#.  -wtT  2  _  -2wt0  0  , 

e  1  (1  -  e  1  )  a  0.5  e  z  ,  u>.  =  — 103 
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giving  T=43  km.  Thus,  the  use  of  commonly  accepted  "realistic"  isostatic  para¬ 
meters  corresponds  quite  closely  to  spherical  harmonic  reference  fields  with 

around  180. 
max 

The  covariance  function  for  the  compensated  model  will  like  (5.6)  be  a  linear 
combination  of  simple  Poisson  covariance  functions 


(5.9) 


C(s) 


‘(77 2^f), 

c  m 


D'=2D,  D^=2DC,  D;=D+DC, 
r  =/s2+D' 2  etc. 


The  potential  covariance  function  K(s)  corresponding  to  (5.9)  will  be  (cf.  Table  1) 


K(s)  *a  (-log(D'+r)  -  log(D^+rc)  +  21og(D(J1+r[n)) 


=  a  log 


(Dm+rm)2 

(D’+r)(D^.+r) 


(5.10) 


which  does  not  exhibit  the  singularities  of  the  individual  "components",  the  geoid 
variance  being 


a2 

n 


K(0)  =  log 


Jmi 

D  D_ 


(5.11) 


The  depths  D 
spectrum  analysis. 


C(0)  =  o2 


Ag 


as 


and  Dc  are  assumed  to  have  been  chosen  e.g.  based  on  a  power 
The  constant  a  are  then  given  from  the  gravity  variance 


C(0) 


(5.12) 


which  by  introduction  of  a  dimensionless  "variance  scaling  factor"  g  yields 


a  =  402C(0)6, 


(5.13) 


The  factor  6  takes  into  account  the  "lost"  power  at  low  frequencies  compared  to 
the  simple  Poisson  model. 
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Returning  to  the  Ohio  example  (Figure  12),  estimates  for  the  geoid  undulations 
associated  with  the  two  source  components  may  now  be  found.  By  inspecting  Figure  12, 
a  fairly  large  amount  of  power  is  seen  to  remain  below  0.5  cycles/degree,  correspon¬ 
ding  to  £=180.  Therefore,  to  take  into  account  also  errors  of  the  reference  model, 
a  fairly  deep  compensation  level  Dc  is  needed.  From  the  position  of  the  3  dB-point 
a  value  of  around  100  km  seems  to  be  reasonable. 

For  Dc=100  km,  the  s-factors  for  the  4.5  km  and  20  km-layers  are  1.01  and 
1.22  respectively.  Therefore  the  23  dB  difference  on  the  DC-values  corresponds 
to  the  slightly  different  variance  ratio  of  0.120,  yielding  r.m.s.  variations 
for  the  shallow  and  deep  components  at  4.0  mgal  and  11.7  mgal .  Using  these  numbers 
and  (5.11)  and  (5.13)  the  r.m.s.  geoid  variations  are  then  found  to  be  5.0  cm 
and  40.5  cm,  yielding  a  total  r.m.s.  geoid  variation  of  40.8  cm.  The  use  of  a 
more  shallow  compensation  level  means  that  the  variance  estimates  contain  less 
of  the  errors  of  the  reference  field.  For  Dc=50  ^m,  r-m*s-  9e°id  variations 
thus  diminish  to  5.0  cm  and  29.5  cm  for  the  shallow  and  deep  sources,  corresponding 
to  a  total  variation  of  29.9  cm.  For  a  crude  comparison  it  may  be  mentioned  that 
9  GPS-determined  geoid  undulations  in  a  smaller  area  of  central  Ohio  (area  ex¬ 
tent-max.  40  km)  have  yielded  a  standard  deviation  of  the  9  points  around  22 
cm  (number  derived  from  undulation  difference  data  given  by  Engelis  et  al . ,  1984.) 

To  complete  the  discussion  of  the  "compensated"  model,  finally  second  order 
derivatives  will  be  considered.  Since  these  quantities  primarily  are  due  to  shallow, 
high-frequency  gravity  field  variations,  little  difference  would  be  expected  between 
the  "simple"  and  "compensated"  Poisson  model.  In  fact,  it  turns  out  that  the  scaling 
factor  8  exactly  opposes  the  compensation  effect,  yielding  identical  variance  ex¬ 
pressions  in  the  two  cases. 

Consider  first  the  vertical  gravity  gradient.  The  "Poisson"  covariance  function 
of  this  quantity  will  be 


G'(s) 


(5.14) 


(cf.  Table  1)-,  and  thus  the  familiar  horizontal  gradient  variance  Gfl  (Moritz  1980, 
p.  177)  is  seen  to  be 


(5.15) 


-.■k"-.'  ■  .'T 
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For  the  compensated  model,  variance  scaling  yields  similarly  by  (5.13) 


G 


o 


=  I  c  ,± 

2  ^0  Q2 


J_n  =  iia 

Q2  2  02 
m 


(5.16) 


i.e.  an  identical  expression. 

In  the  Ohio  example,  the  shallow  and  deep  sources  thus  correspond  to  r.m.s. 
gradient  variations  (/G^)  of  10.9  E  and  7.1  E  respectively  (1  E=10-4mga1/ m) , 
with  a  total  r.m.s.  variation  of  13.0  E.  This  number  may  again  be  compared  with 
observation  data.  Ca.  300  torsion  balance  stations  exist  in  an  area  of  south¬ 
western  Ohio  (surveyed  by  Badekas  of  O.S.U.  nearly  two  decades  ago),  and  for  these 
stations  a  /CT”  -value  of  c.  18  E  was  found  (Tscherning,  1976,  Figure  2a).  The 
variances  compare  reasonably  well,  especially  considering  the  limited  resolution 
of  the  gravity  data  (4‘  grid)  underlying  the  power  spectrum  plot  of  Figure  12. 

The  combined  results  of  the  Ohio  example  used  in  this  section  are  shown 
in  Table  2.  Deflections  of  the  vertical  follow  simply  from  the  gravity  results 
using  the  conversion  factor  6.7  mgal/arcsec  of  section  2.1.  A  more  thorough 
3-layer  analysis  of  part  of  this  example  area  will  be  included  in  the  next  section. 


Table  2 

Two-Component  Covariance  Model  for  Ohio  4°x4°area. 

(white  noise  compensated  density  layer  representation,  Rapp-180  reference  field) 


Source 

Depth 

(km) 

C(0) 

(mgal 2) 

r.m.s.  variation  of 

geoid 

(m) 

deflections 

(arcsec) 

free-air 
gravi ty 
(mgal ) 

gravity 

gradient 

(E) 

"sedimentary" 

4.5 

16.3 

0.05 

0.6 

4.0 

10.9 

"deep  crustal" 

20 

136.7 

0.40 

1.7 

11.7 

7.1 

total 

— 

153.0 

0.41 

1.8 

12.4 

13.0 

_  a  *»  •  ,  a  a  .  .  _  «  ’  »  *  »  '  •  “J *  *  •  *  “  •"  *"  .*  •*-  **  «*.  ^  '  a'  «  *  • 
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6.  Examples  of  Local  Empirical  Covariance  Functions  and  Power  Spectra 

In  this  section  some  examples  will  be  given  on  the  shape  of  local  power 
spectra  and  covariance  functions  for  gravity  anomalies  in  selected  areas  of  the 
United  States.  Areas  of  fairly  dense  anomaly  coverage  have  been  chosen,  in  order 
to  get  the  best  estimates  of  the  power  spectra  so  that  possible  "straight-line 
segments"  associated  with  major  geologic  interfaces  would  show  up  more  clearly. 
Also,  areas  of  dense  data  coverage  permits  the  degree-variances  to  be  derived 
with  confidence  up  to  high  degrees,  in  the  present  study  to  t=5400. 

The  data  used  -  free-air  anomalies  and  terrain-corrected  Bouguer  anomalies  - 
were  supplied  by  the  National  Geodetic  Survey.  In  order  to  work  with  local  co- 
variance  functions,  the  "Rapp-180"  spherical  harmonic  expansion  of  the  geopoten¬ 
tial,  complete  to  degree  and  order  180  (Rapp,  1981),  has  been  subtracted  for 
all  data  to  yield  residual  anomalies.  Furthermore,  since  one  of  the  primary 
interests  is  in  possible  "geologic"  effects  in  the  power  spectrum,  the  influence 
of  the  topography  has  been  removed  in  the  mountainous  areas.  The  elimination 
of  the  topographic  effects  on  the  gravity  anomalies  has  been  done  using  a  re¬ 
sidual  terrain  model  (RTM)  reduction,  where  topographic  irregularities  relative 
to  a  smooth  mean  elevation  surface  are  computationally  removed.  As  mean  eleva¬ 
tion  surface  a  spherical  harmonic  expansion  of  the  topography,  corresponding 
to  the  gravity  expansion,  has  been  used  (Rapp,  1982). 

The  RTM-reduction  is  very  convenient,  since  it  effectively  corresponds  to 
a  removal  of  all  topographic  effects  above  t=180,  while  nothing  in  principle 
is  done  below  this  degree  (the  "Rapp-180"  gravity  field  also  contains  effects 
of  the  topography).  As  shown  in  (Forsberg,  1984),  the  RTM  reduction  corresponds 
usually  within  a  fraction  of  a  mgal  to  a  Bouguer  reduction  to  the  reference  ele¬ 
vation  level.  Thus,  the  residual,  terrain-reduced  gravity  anomalies  used  here 
may  be  written 


4s  ■  AgBA  -  ag8A  *  AgBA  -  (ag18t  -  a-**,,,) 


(6.1) 


with 


*9BA 

A9180 


RTM  anomaly,  corresponds  to  Bouguer  anomalies  above  t*180 
Conventional,  terrain-corrected  Bouguer  anomaly 
Free-air  anomaly  of  the  "Rapp-180"  field 
Elevation  from  the  180  x  180  topography  expansion 


For  the  practical  computations,  Ag18Q  and  h18Q  were  computed  in  a  0.25°x0.25° 
grid,  from  which  values  at  individual  stations  were  determined  by  a  simple  bi¬ 
linear  interpolation  scheme. 

Gravity  anomalies  for  four  different  areas  have  been  investigated.  These 
areas  will  in  the  sequel  be  designated  by  the  name  of  the  state  inside  which 
they  are  located: 

"OHIO":  2°x2°  lowland/moderately  hilly  area  SW  of  Columbus, 

1  a t 38°-40°  N,  Ion.  85°-83° W  (not  identical  to  the  Ohio 
example  of  the  last  section). 

"COLORADO":  2°x  2°  alpine  area  of  the  Rocky  Mountains,  W  of  Denver 
and  Colorado  Springs,  lat.  38°-40°N,  Ion.  107°-105°W. 

“CALIFORNIA":  2°x  2°  mixed  area,  containing  the  central  parts  of  the 

Sierra  Nevada  mountains  and  a  part  of  the  California  valley, 
lat.  36° -38°  N,  Ion.  120°-118°N. 

"NEW  MEXICO":  4°x  4°  mountainous  area,  covering  most  of  New  Mexico, 
including  the  "White  Sands"  area.  lat.  32°-36°N,  Ion. 

109°-105°  W. 

Except  for  New  Mexico,  these  areas  represent  new,  more  detailed  investigations 
than  the  results  previously  presented  (in  a  somewhat  different  context)  in  Forsberg 
(1984).  The  2°x  2°  areas  have  been  chosen  to  secure  a  sufficiently  dense  coverage 
of  gravity  anomalies,  to  allow  a  fairly  dense  grid  (2‘ x  2.5' )  of  gravity  anomalies 
to  be  predicted.  An  example  of  the  data  coverage  is  shown  in  Figure  15. 

The  gridding  procedure  and  analysis  were  similar  to  the  method  of  Forsberg 
(1984).  First,  gravity  anomalies  were  screened,  keeping  only  one  anomaly  per 
1  * x  1*  "pixel".  The  screened  anomalies  were  then  gridded  to  a  2‘x2.5'  (-3.7x3. 6  km) 
grid  using  a  truncated  collocation  algorithm,  where  a  grid  point  value  is  predicted 
from  the  five  closest  observations  only.  The  dense  prediction  grid  necessitates 
some  caution  when  interpreting  results  for  higher  wave  numbers:  Although  the 
data  coverage  for  the  areas  is  fairly  dense,  it  is  evident  from  Figure  15  that 
the  data  might  be  too  sparse  in  parts  of  the  areas,  yielding  some  influence  of 
the  grid  prediction  method  on  the  shape  of  the  power  spectrum  at  high  wave  numbers. 

From  the  gridded  RTM  gravity  anomalies  power  spectra,  covariance  functions 
and  degree  variances  were  constructed  using  the  two-dimensional  Fast  Fourier 
Transform  (FFT).  The  power  spectrum  is  obtained  from  the  Fourier  transform  by 
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Figure  15  Gravity  data  coverage 
in  the  2°x2°  area  of  California. 
The  California  valley  is  to  the 
lower  left,  while  the  Sierra  Nevada 
mountains  run  approximately  NVI-SE 
through  the  block.  Totally  3172 
1‘xl'  pixels  contain  gravity 
data. 


and  by  an  inverse  transform  the  two-dimensional  covariance  function  is  obtained 


C(x,  y)  «  Fourier-1  Ugg(u,  v)) 


(6.3) 


To  yield  "isotropic"  quantities,  $  and  C  were  subsequently  averaged  along  circles 
of  constant  w  */uT+v7'and  s*/xz+yz  respectively,  and  from  the  averaged  power 
spectrum  potential  degree  variances  were  finally  derived  by  (2.9)  and  (2.34): 


R2 


U-l)2 


1  l+h 
2tt  (i-l)2 


(6.4) 


The  spectrum  &g  estimated  using  the  fast  Fourier  transform  will  be  affected 
by  errors  due  to  the  finite  extent  of  the  area  in  question,  and  some  suitable 
window  filter  should  ordinarily  be  used.  In  the  present  investigation,  however. 


windowing  turned  out  to  be  a  disadvantage  due  to  the  unavoidable  loss  of  power 
in  the  spectrum  -  and  through  the  radial  averaging  procedure  already  some  smoothing 
is  done.  Window  filters  have  therefore  not  been  applied  here,  resulting  in  somewhat 
"wiggling"  spectra. 

Results  for  the  2°x  2°  areas  are  shown  in  Figures  16-18.  Each  figure  shows 
from  top  to  bottom  graphs  of  o^,  ^(w)  and  C( s )  shown  with  a  thick,  solid  line. 

The  degree  variances  o  are  normalized  and  shown  logarithmic,  with  the  global 
degree- variance  models  of  Kaula  (3.35)  and  Tscherning-Rapp  (3.36)  also  shown 
for  comparison  purposes.  The  power  spectrum  is  shown  logarithmic  in 
dB  (dB=10  lO9104gg)»  with  unit  (0  dB)  1  mga12degree2.  Wave  numbers  are  shown 
in  cycles/degree  (-^).  with  maximum  15  cycles/degree  being  the  Nyquist  frequency 
for  a  2'  data  grid.  With  broken  lines  are  shown  arbitrary  tentative  straight-line 
segments,  with  depth  D  determined  from  the  slope  of  the  line.  Finally,  the  covariance 
function  C(s)  is  shown  in  units  of  mgal2.  For  reference  purposes,  the  covariance 
function  has  been  fitted  with  Moritz'  two  "basic"  planar  covariance  functions 
(3.36)  and  (3.27):  the  reciprocal  distance  covariance  (R)  and  the  Poisson  covariance 
(P),  shown  with  thin  line  also  in  the  power  spectrum  plot.  The  covariance  functions 
were  fitted  simply  by  requiring  variance  C0  and  correlation  length  x,  to  match 

d. 

the  empirical  values.  To  illustrate  possible  anisotropy,  the  central  part  of 
the  two-dimensional  covariance  function  C(x,  y)  is  shown  as  a  contour  plot  to 
the  right,  with  a  contour  interval  of  0.2  C0.  For  a  perfectly  isotropic  process 
this  plot  would  show  a  series  of  concentric  circles. 

(The  corresponding  plot  of  the  New  Mexico  area,  derived  from  a  more  coarse 
4 ' x  5 '  anomaly  grid,  may  be  found  in  (Forsberg,  1984,  Figure  34)). 

Inspection  of  the  power  spectra  shows  that  straight-line  segments  are  not 
very  conspicius.  Indeed  the  spectra  seems  to  be  of  more  or  less  the  same  overall 
shape.  The  depths  assigned  to  various  parts  of  the  power  spectra  in  the  plots 
must  be  viewed  with  great  caution,  since  they  are  characterized  by  very  large 
uncertainties.  Especially  for  the  higher  wave  numbers  error  sources  due  to  the 
gridding,  aliasing,  and  errors  in  the  anomalies  and  terrain  reductions  play  quite 
a  substantial  role  in  shaping  the  spectrum.  The  influence  of  random,  uncorrelated 
errors  in  the  grid  points  may  be  ascertained  fairly  easily  using  Parseval's  equa¬ 
tion  (2.16):  Viewing  the  errors  as  bandlimited  white  noise  (in  a  frequency  square 
of  side  length  2x15  cycles/degree),  it  follows  that  an  error  variance  of  1  mgal2 
corresponds  to  -30  dB  in  Figures  16-18.  Thus  e.g.  the  tail  "d-  1  km"  of  the 
Colorado  spectrum  (Figure  17)  may  be  explained  alternatively  as  the  effect  of 
uncorrelated  gravity  grid  errors  of  variance  16  mgal2  (-  -18  dB  in  the  plot),  not 
too  unrealistic  considering  the  extreme  alpine  topography  of  the  area. 
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Figure  16  Degree  vari¬ 
ances,  power  spectrum 
and  covariance  function 
for  free-air  gravity 
anomalies  in  a  2°x  2° 
area  of  Ohio.  See  text 
for  detailed  figure 
explanation. 
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Figure  18  Degree-vari- 
ances,  power  spectrum 
and  covariance  function 
for  residual  Bouguer 
(RTM)  anomalies  in  a 
2°x  2°  area  of  California 
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The  geophysical  significance  of  the  Indicated  depth  values  is  probably  limited, 
the  depth  values  possibly  representing  merely  tendencies  to  enhanced  density  vari¬ 
ations  at  certain  depth  ranges  rather  than  undulations  of  some  density  contrast 
interface.  The  lack  of  clear  "layering"  in  the  power  spectrum  is  as  earlier  men¬ 
tioned  typical  of  gravity  problems,  opposed  to  many  magnetic  examples.  The  sources 
of  the  gravity  field  variations  (besides  the  topography)  tends  to  be  distributed 
throughout  the  crust  and  upper  mantle,  in  a  random  fashion  quite  similar  for  many 
different  types  of  geologic  settings. 

Irrespectively  of  the  geophysical  interpretation  of  the  indicated  depths, 
they  may  be  used  to  provide  improved  multi-layer  local  covariance  functions  of 
the  gravity  field  in  the  respective  areas.  For  each  of  the  four  different  areas 
three  depth  values  have  been  selected  as  outlined  in  the  figures,  and  a  total 
covariance  function  is  built  up  by  using  compensated  Poisson  models  at  each  of 
the  specified  depths.  This  gives  three  unknowns  -  the  variance  associated  with 
each  layer  -  but  since  the  total  variance  is  given  from  the  data,  only  two  free 
parameters  are  left.  In  addition  also  the  compensation  depth  Dc  may  be  varied, 
especially  in  order  to  model  the  apparent  power  left  below  degree  180,  originating 
in  deficiencies  in  the  used  spherical  harmonic  expansions  of  the  gravity  field 
and  the  topography. 

The  covariance  models  have  in  the  present  study  been  fitted  to  the  power 
spectra  using  a  simple  interactive  ad  hoc  technique,  utilizing  a  graphical  terminal. 
Each  fitted  power  spectrum  is  of  form  (5.6) 

*gg(“)  =  27rIx  Ve~2“0i  +  e’2wDc  '  2e'“(Di+Dc))  (6.5) 


where  is  transformed  to  a  gravity  component  variance  by  (5.12).  Given  the 
variances  the  geoid  undulation  variance  and  gravity  gradient  variance  are  sub¬ 
sequently  found  by  (5.11)  and  (5.16).  The  results  are  shown  in  Table  3  and  Figure  19 
The  covariance  functions  of  Table  3  are  examples  of  "tailored"  covariance 
functions,  applicable  e.g.  for  optimal  error  studies,  simulations  etc.  Since 
they  have  a  good  fit  in  the  entire  spectral  range,  they  will  be  good  both  for  geoid 
and  second  order  derivative  applications.  The  models  show  clearly  the  general 
feeling  of  gravity  field  variations:  the  geoid  is  sensitive  to  deep  density  vari¬ 
ations,  while  second-order  gradients  primarily  reflect  shallow  sources.  Comparing 
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Flgure  19  Examples  of  3-component  compensated  Poisson  power  spectra  models  for 
the  sample  areas.  Derived  by  a  simple  graphical  fit,  constrained  by 
conservation  of  the  total  variance.  Parameters  listed  in  Table  3. 

the  2°x  2°  Ohio  area  to  the  previous  example  of  Figure  12  (4°x4°  area),  it  Is 
seen  that  quite  a  substantial  part  of  the  gradient  variations  was  not  resolved 
In  the  example  of  the  last  section.  The  Increased  resolution  in  the  gravity  grid 
of  the  present  example  thus  provides  for  a  better  agreement  with  the  torsion  balance 
results  mentioned  in  the  last  section. 


Table  3 


3-Layer  Compensated  White  Noise  Density  Variance  Models 
for  Residual  Gravity  In  Sample  Areas  of  the  United  States 
(Topography  Removed  by  RTM- Reduct ion.  Except  for  Ohio) 


Area 

Location 

Layer 

Depth 

(km) 

Gravity 
Variance 
(mgal 2 ) 

r.m.s.  variation  of 

geo  id 

undulations 

(m) 

gravity 
anomalies 
(mgal ) 

horizontal 
grav.  gradient 
(E) 

— 

2 

■RQI 

3.6 

22.3 

OHIO 

4.5 

4.3 

11.7 

20 

210.3 

14.5 

8.9 

(Dc=50  km) 

total 

242.1 

0.38 

15.6 

26.7 

1 

26.1 

0.02 

5.1 

62.6 

COLORADO 

3.3 

24.0 

0.05 

■SB 

18.2 

16 

322.8 

0.54 

mm 

13.8 

(0C*100  km) 

total 

373.0 

0.54 

19.3 

66.6 

2 

30.8 

0.04 

5.5 

34.0 

CALIFORNIA 

4 

46.4 

0.08 

6.8 

20.9 

26 

1148.8 

1.37 

33.9 

16.0 

(Dc=100  km) 

total 

1226.0 

1.37 

35.0 

42.9 

2.5 

22.9 

0.04 

4.8 

23.4 

NEW  MEXICO 

12 

62.6 

0.20 

■fifl 

8.1 

40 

112.4 

0.55 

■a 

3.2 

(Dcs100  km) 

total 

197.8 

0.58 

14.1 

25.0 

With  the  apparent  quite  high  similarity  between  the  various  power  spectra, 
it  Is  naturally  then  to  ask  the  question:  What  is  the  "best"  overall  shape  of 
the  covariance  functions  in  terms  of  the  "simple"  covariance  models?  As  earlier 
stated,  this  is  in  general  believed  to  be  the  logarithmic  clays  of  the  covariance 
functions.  That  this  is  indeed  the  case  may  be  seen  from  Figure  20,  showing  the 
derived  degree  variances  on  a  double  logarithmic  plot,  power  functions  a“2,  a'3 
and  t”4  have  additionally  been  plotted  for  reference  purposes.  Remembering  the 
asymptotic  forms  of  some  of  the  "elementary"  covariance  functions  of  Table  1. 

poisson  :  -  a-1 

reciprocal  distance:  a  -  a-2 

X> 

logarithmic  :  -  a-3 

it  is  dear  that  the  logarithmic  types  should  be  preferred  in  general.  The  simple 
planar  logarithmic  covariance  function  of  Table  1  is  singular,  so  modifications 
analagous  to  the  compensated  Poisson  model  will  be  needed  to  eliminate  the  influence 
of  the  low  wave  numbers.  Alternatively,  the  well-known  spherical  Tscherning-Rapp 
model  (3.36)  may  be  used,  with  adaption  to  local  applications  obtained  through 
a  change  in  variance  factor,  Bjerhammar  sphere  depth  and  removal  of  lower  degree- 
variances  corresponding  to  the  spherical  harmonic  reference  field. 

The  degree  variances  of  Figure  20  are  seen  to  be  systematically  too  low  '*in- 
pared  to  Kaula's  rule  for  the  higher  degrees.  Kaula's  rule  represents  a  kind 
of  global  average  of  the  actual  gravity  field  variation,  including  the  effects 
of  the  topography. 

For  a  valid  comparison,  the  topographic  effects  must  therefore  not  be  removed 
from  the  basic  gravity  grid  used.  In  Figure  21  are  shown  the  corresponding  em¬ 
pirical  degree  variances  for  unreduced  gravity  data,  i.e.  free-air  anomalies, 
derived  from  the  gridded  Bouguer  anomalies  by  an  inverse  Bouguer  reduction  using 
elevations  from  detailed  ( 0 . 5 ' xO . 5 ' )  digital  terrain  models  of  the  areas.  Now 
the  spread  in  the  degree  variances  is  much  larger,  reflecting  the  differences 
in  terrain.  The  logarithmic  asymptotic  behaviour  is  still  clear,  and  now  Kaula's 
rule  performs  somewhat  better,  indicating  e.g.  that  New  Mexico  is  a  quite  typical 
area  on  a  global  basis. 

The  logarithmic  nature  of  empirical  gravity  field  covariance  functions  may 
be  used  for  some  general  implications  on  the  possible  types  of  density  distributions. 
If  we  e.g.  accept  the  stationary  thin-layered  earth  model,  where  the  earth  was 
described  by  a  number  of  thin,  independent  layers,  the  total  power  spectrum  was 
found  in  section  (3.2)  to  be  of  form 


POT-  DEGREE  VARIANCE  (LOGIO) 
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,U)  =  l  (Z^f  ^p(u>)  e-2“di 


(6.6) 


where  az.  is  the  layer  thickness,  d,  the  depth  and  (j*1  the  density  power  spectrum 
i  i  pp 

of  the  layer.  The  form  of  the  power  spectrum  for  infinitely  thin  layers  will  thus  be 


>n>)  -  f  z)e"2w2dz 

99  o  pp 


(6.7) 


If  the  layers  are  assumed  to  be  identical ,  a  logarithmic  power  spectrum  (i.e., 
4>gg(w)  -  u>~2)  requires 


tp  (w)  -  — 
pp  <*> 


(6.8) 


i.e.  a  red  noise  distribution,  where  density  anomalies  of  large  spatial  extent 
are  on  the  average  of  bigger  magnitude  than  the  density  features  of  less  extent, 
a  quite  reasonable  statistical  model  rather  than  the  "strong"  white  noise  assumption. 

If,  on  the  other  hand,  the  concept  of  white-noise  density  layers  is  maintained, 
it  then  follows  that  the  density  variance  a2  must  be  given  by  an  inverse  Laplace 

D 

transform 


►flflU)  -  4  -  r  o2(z)e‘2wZdz 

gg  o  p 


(6.9) 


yielding 


»p<z>  -  z 


(6.10) 


i.e.,  the  density  anomaly  variance  increases  linearly  with  depth,  a  very  unphysical 
result  considering  our  general  ideas  of  the  earth's  interior. 
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7.  Stannary  and  Conclusions 

In  the  present  report  various  topics  relating  to  covariance  functions  and 
density  distributions  have  been  treated.  Emphasis  has  been  on  local  covariance 
functions,  where  a  high  degree  and  order  spherical  harmonic  reference  field  are 
substracted  prior  to  modelling  and  analysis.  By  working  with  the  residual  field, 
a  number  of  simplifying  assumptions  may  be  used,  such  as  the  flat  earth  formulation 
and  the  identification  of  gravity  anomalies  with  gravity  disturbances.  With 
these  assumptions  it  was  e.g.  shown  how  collocation  may  be  interpreted  in  a  very 
simple  fashion  in  terms  of  generalized  point  mass  modelling. 

Apart  from  the  collocation  interpetation ,  the  various  topics  of  the  present 
report  may  be  generalized  easily  using  a  spherical  earth  formulation.  An  example 
(which  actually  was  part  of  the  inspiration  for  the  work  of  the  present  report) 
may  be  found  in  Sunkel  (1981),  who  considers  the  white  noise  density  layer  covari¬ 
ances  on  a  strictly  spherical  basis.  For  local  applications,  however,  spherical 
results  follow  much  more  easily  by  the  simple  relationship  existing  between  the 
planar  power  spectrum  and  the  spherical  harmonic  degree  variances. 

This  relationship  was  elaborated  upon  in  Section  2,  in  addition  to  the  most 
basic  description  needed  for  using  the  planar  formulation. 

In  Section  3  then  the  relationship  between  the  well  known  covariance  functions 
(Poisson,  Logarithmic  etc.)  and  "mass  bodies"  characterized  by  white-noise  distribu¬ 
tions  of  density  anomalies,  dipole  moments  etc.,  were  derived  using  the  ensemble 
averaging  theorem,  well  known  to  geophysicists  as  the  method  of  Spector  and  Grant 
(1970).  The  use  of  the  theorem  opens  in  principle  for  a  direct  link  between 
the  shape  of  the  power  spectrum  (and  thus  the  covariance  function)  and  depth 
to  "disturbing"  interfaces  in  the  earth's  interior. 

The  simple  covariance  functions,  described  e.g.  as  generated  by  a  white-noise 
density  layer,  vertical  mass  lines  etc.,  turned  out  to  be  themselves  expressable 
as  "gravity”  effects  of  similar  (but  not  the  same)  bodies,  only  at  the  double 
depth.  This  was  used  for  the  side  remark  in  Section  4,  the  collocation  interpreta¬ 
tion.  It  was  e.g.  found  simply  that  collocation,  using  the  Poisson  covariance 
function  corresponds  to  a  point  mass  modelling,  where  a  point  mass  should  be 
located  below  each  gravity  observation, a  vertical  mass  line  below  each  geoid 
observation  point  etc.,  at  depths  two  times  the  depth  to  the  Bjerhammar  sphere. 

For  practical  local  covariance  functions,  in  accordance  with  geologically 
reasonable  density  distributions,  multi-layer  models  of  independent  density  vari¬ 
ations  will  be  necessary,  corresponding  to  different  Bjerhammar  sphere  depths 
for  the  different  constituents.  In  Section  5  especially  the  multi-layer  Poisson 


vs 


model  was  treated  In  some  detail,  and  a  modification  based  on  a  pseudo-isostatically 
"compensated  white  noise  density"  model  was  proposed  to  model  the  low  frequency 
behaviour  of  the  residual  gravity  field  more  properly.  It  turned  out  that  with 
a  perfect  180x  180  spherical  harmonic  reference  field,  the  "optimum"  compensation 
level  was  only  slightly  deeper  than  the  commonly  accepted  compensation  depths 
of  Airy  isostasy.  The  expressions  for  the  covariances  etc.  of  the  compensated 
model  turned  out  to  be  very  simple,  since  the  formulas  were  obtained  as  a  simple 
linear  combination  of  the  original  Poisson  formulas. 

Finally,  in  Section  6,  actual  power  spectra  for  the  various  sample  areas 
were  analyzed,  with  the  primary  object  to  detect  the  existence  of  possible  dis¬ 
turbing  interfaces.  The  geophysical  significance  of  the  results  turned  out  to  be 
questionable,  an  experience  obtained  by  many  geophysicists  having  tried  to  use 
statistical  interpretation  techniques  for  gravity  anomalies.  However,  the  in¬ 
ferred  layer  parameters  are  very  useful  in  providing  optimal  "tailored"  covari¬ 
ance  models,  and  numerical  values  for  such  models  were  given  for  the  different 
sample  areas,  and  illustrated  by  providing  statistics  for  the  geoid  and  second 
order  gradient  variances. 

Forgetting  the  multi-layer  models,  the  overall  shape  of  the  covariance  func¬ 
tion  turned  out  to  be  logarithmic,  with  degree  variances  decaying  like  jT3  up 
to  at  least  degree  5400.  Kaula's  rule  turned  out  to  be  quite  realistic  for  the 
"actual"  gravity  field,  while  the  power  was  significantly  too  high  when  terrain- 
reduced  data  were  used. 

Although  the  statistical  "covariance  interpretation"  techniques  have  not 
too  bright  prospects,  I  think  it  could  be  of  interest  for  future  studies  to  take 
more  sample  areas,  especially  areas  in  geologic  settings  where  a  "white-noise 
layer"  description  would  be  expected  to  have  some  validity,  e.g.  for  describing 
undulations  in  salt/sediment  Interfaces  In  salt  dome  provinces  (like  e.g.  the 
Gulf  coast  of  Texas  or  northern  Jutland,  Denmark).  In  such  cases  important  in¬ 
formation  on  the  very  local  variability  and  stationarity  of  the  "detailed"  co- 
variance  functions  could  additionally  be  derived.  A  prerequisite  for  such  studies 
would  be  a  very  dense  gravity  coverage,  combined  with  power  spectrum  estimation 
techniques  yielding  reliable  estimates  even  for  small  sample  areas,  inside  which 
the  "layering"  assumption  is  reasonably  valid. 
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